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Abstract 

F[igh index differential algebraic equations (DAEs) are ordinary differential equations (ODEs) 
with constraints and arise frequently from many mathematical models of physical phenomenons 
and engineering fields. In this paper, we generalize the idea of differential elimination with Dixon 
resultant to polynomially nonlinear DAEs. We propose a new algorithm for index reduction of 
DAEs and establish the notion of differential algebraic elimination, which can provide the differ¬ 
ential algebraic resultant of the enlarged system of original equations. To make use of structure of 
DAEs, variable pencil technique is given to determine the termination of differentiation. More¬ 
over, we also provide a heuristics method for removing the extraneous factors from differential 
algebraic resultant. The experimentation shows that the proposed algorithm outperforms existing 
ones for many examples taken from the literature. 

Key words: Index reduction; Differential algebraic resultant; Variable pencil; Differential 
algebraic equations 


1. Introduction 


Modeling with differential algebraic equations (DAEs) plays a vital role in a variety of ap¬ 
plications Olhll . for constrained mechanical systems, control theory, electrical circuits and chem¬ 
ical reaction kinetics, and many other areas. In general, it is directly numerical computations 
difficult to solve the system of DAEs. The index of DAEs is a measure of the number of times 
needed to differentiate it to get its equivalent low index or ordinary differential equations (ODEs). 
There exist many different index concepts for the specific DAEs, such as the differentiation in¬ 
dex iSli, perturbation index §[11, tractability index |[H . structural index ll^ . and Kro- 
neckerindex 
DAEs ifliilQ 


l3in . There has been considerable research for the general linear and low index 
21L 31 [. In particular, it may only solve some special DAEs by the directly nu¬ 
merical methods |l Mll8n . It is more difficult to solve the system of high index nonlinear DAEs 


ijiliSlilllll El . Therefore, index reduction techniques may be necessary to get a solution 

m- 
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Index reduction in the pre-analysis of DAEs solving is an active technique of research. It 
is equivalent to applying a sequence of differentiations and eliminations to an input system of 
DAEs. In ll^ . Pantelides gave a systematic way to reduce the high index DAEs to low index 
one, by selectively adding differentiated forms of the equations already appear in the system. 
However, the algorithm can succeed yet not correctly in some instances and be just first order 
1 2711 . Campbell’s derivative array theory needs to be computationally expensive especially for 
computing the singular value decomposition of the Jacobian of the derivative array equations 
using nonlinear singular least squares methods il. In ii, Mattsson et al. proposed the dummy 
derivative method based on Pantelides’ algorithm for index reduction, which is an algebraic view¬ 
point. In ll23l] . Pryce proposed the signature matrix method (also called E-method), which can 
be viewed as an extension of Pantelides’ method for any order. Recently, Wu et al. generalized 
the E-method for DAEs to partial differential algebraic equations with constraints (PDAEs) 

Qin et al. presented the structural analysis of high index DAEs for process simulation by the 
E-method 12411 . But the E-method relies heavily on square (i.e. the same number of DAEs and 
dependent variables) and sparsity structure, which is confronted with the same drawback that can 
succeed yet not correctly in some DAEs arising from the specific types. 

A principal aim of this paper is the development of an efhcient differential elimination ap¬ 
proach for index reduction of DAEs that extends the direct elimination treatment of S. Erom 
the algebraic standpoint, differential elimination algorithms which are key for simplifying sys¬ 
tems of polynomially differential equations and computing formal power series solutions for 
them. The underlying theory is the differential algebra of Ritt |28] and Kolchin ifisll . Differen¬ 
tial elimination algorithm in algebraic elimination theory is an active field and powerful tools 
with many important applications [ 33 ,[lit]. Almost all of the authors focus on the 

differential elimination theory for ODEs. Only Reid et al. presented an effective algorithm for 
computing the index of polynomially nonlinear DAE and a framework for the algorithmic anal¬ 
ysis of perturbed system of PDAEs. This underlies the jet space approach based on differential 
geometry. 

In this paper, we want to promote the efficient differential elimination algorithm as natural 
generalization of DAEs, which is a direct and elementary approach. In particular, differential 
elimination with Dixon resultant can be solved by eliminating serval variables at a time, sim¬ 
plifying the system with respect to its constraints, or determining its singular cases 13 ■ We 
can directly transform the system of DAEs to its equivalent ODEs by differential algebraic elim¬ 
ination. Differential algebraic elimination is to apply a finite number of differentiations and 
eliminations to uncover all hidden constraints of system of DAEs. We define the new minimum 
differentiation time, which is the weak differentiation index for DAEs/ODEs. It can be used as a 
unified formulation of differentiation times for differential elimination of ODEs and differential 
algebraic elimination of DAEs. Meanwhile, we provide the new index reduction with variable 
pencil and the notion of differential algebraic resultant. In order to overcome the drawback of 
factoring a large polynomial system i34ll . we consider a heuristics method for removing the ex¬ 
traneous factors from the differential algebraic elimination matrix. Our algorithm is also suitable 
for the non-square nonlinear DAEs/ODEs. To the best of our knowledge, it is the first time that 
the generalized Dixon resultant formulation has been directly extended to the system of DAEs. 

The rest of the paper is organized as follows. Section 2 gives a brief description of the gen¬ 
eralized Dixon resultant formulation, and analyzes the size of Dixon matrix and the complexity 
of computing the entries of Dixon matrix. Section 3 proposes the new index reduction proce¬ 
dure for the system of DAEs and defines the weak differentiation index. Section 4 provides the 
differential algebraic elimination algorithm and some basic properties of differential algebraic 
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resultant. Section 5 presents some specific examples in detail and comparisons of our algorithm 
for the system of ODEs. The final section concludes this paper. 


2. Generalized Dixon resultant formulation 

Following Kapur et al. iSil we introduce the concept of generalized Dixon 

resultant formulation and its properties. This technique will play a central role in our subsequent 
analysis. Let X - {xi, X 2 , ■ ■ ■ , x,,] and X = {xi, X 2 , ■ • • , Xn) be two sets of n variables, respectively. 
The determinant of a square matrix A is denoted by det(A). 


Definition 2.1. Let T — {/i,/ 2 , • • ■ ,fn+i} c Q[X] be a set of n + 1 polynomials in n variables. 
The cancellation matrix Cf off is the (n + 1) X (n + 1) matrix as follows: 


Cr 


Mxi,X2,--- ,x„) 
Mxi,X2,--- ,X„) 

fl{Xl,X2,--- ,Xn) 


fn+\{Xi,X2,--- ,X„) 
fn+\ (-^ 1 1 X2^ ■ ' * 5 xf) 
fn+\{X\,X2,--- ,X„) 


/l(Xi,X2, ,X„) 


fn+\{Xi,X2,--- ,X„). 


where f{x\,X2, ■ ■ ■ , x^+i, x<:+ 2 , • ■ • , x„) stands for uniformly replacing xj by xjfor all \ < j < 

k < n in f. The Dixon polynomial off is denoted by Of £ QK 


det{Cf) 


( 1 ) 


the row vector of Dixon derived polynomials off is denoted by Pf, and Dixon matrix off is 
denoted by Df as follows. 


Of = PfVgiOf) = Vx(9T)DrVx(0r), 


( 2 ) 


where YgiOf) is a column vector of all monomials in X which appears in Of, and VxiOf) is a row 
vector of all monomials in X which appears in Of. The determinant of Df is called the Dixon 
resultant, denoted by res{f \, / 2 , ■ • • , fi+f)- 

It is well known that Dixon resultant is a projection operator whose vanishing is a necessary 
condition for the system f to have a common affine solution. However, the Dixon matrix may 
be non-square then its determinant cannot be directly computed. Even if it is square, the Dixon 
resultant vanishes identically without providing any information for the affine solutions. In S, 
Kapur et al. presented a heuristic method to remedy the drawback by extracting a non-trivial 
projection operator. 

Lemma 2.1. ( U^l ) If there exists a column which is linearly independent of all other columns 
in Df, then the determinant of any non-singular rank submatrix ofDf is a non-trivial projection 
operator. 


Remark 2.2. From Lemma 12.71 this method may fail if there is no column which is linearly 
independent of all other columns in Df. However, the method is quite efficient and practical 
as demonstrated in [l4 M 36 HI/, and such failure is very rare even never occurred on the 
numerous problems. Furthermore, the projection operator may contain extraneous factors in the 
Dixon resultant. 
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In this article, we shall use the following properties of Dixon resultant. 

Lemma 2.3. Dixon resultant can be expressed as a linear combination of original poly¬ 

nomial system ‘T, 

n-v\ 

res{f\ , /2, • • • , /„+i) = 'Yj 

i=l 

where Ki is a polynomial with respect to X and can be deduced from Pjr. Moreover, it has been 
proved that the extraneous factors mentioned above may include three parts which are taken from 
Pyr, Dyr and the resulting resultant expression by substituting Pjr, respectively. 

Remark 2.4. Extraneous factors arising from Dixon resultant is a troublesome problem when it 
is used for elimination in a variety of applications. Gather-and-Sift method 4331/ is a complete 
method to remove extraneous factors by the simplicial decomposition algorithm. But it suffers 
from very high computational complexity because of the intermediate expression swell in sym¬ 
bolic computation. Therefore, we mainly use the technique based on Lemma 12.31 which can be 
viewed as a heuristic method. 


Theorem 2.5. The size of Dixon matrix is at most n\ H/Li d-i x n\ HILi di> the complexity 
of computing the entries of Dixon matrix is 0(d^(n\ Yi’Lidi)^) in worst case, where dj is the 
highest degree of variable x,-. 


Proof. Similar to the proof of computing the entries of Dixon matrix in the bivariate case and 
combine with the multivariate Sylvester resultant and the general case in 03611381] . □ 


Remark 2.6. Here we give the size and computational complexity of Dixon matrix in the general 
setting. In particular, the complexity of computing the entries of Dixon matrix is a new result. 
The highest degree di of variable x,- can be obtained by using the algorithm in i2^/ . To make use 
of sparsity in polynomial systems, bound on the size of Dixon matrix of the specific systems is 
derived in terms of their Newton polytopes in Ml- 


3. Index reduction algorithm 

In this section, let 6 - djdt denote the differentiation operator, let be a differential ring, i. 
e., a commutative ring with unit and a differentiation 5 acting on it. Let No = {0,1, • • ■ 

1 / 6 Nq represents a multi-index v - {vi,V2, - ■■ ,v„ )^meN,Y ^ {yuy2,.-.,y„,}- If r 6 No, 
then order of 5'' is ord{6'') - r, we denote y® the A:-th derivative of yj and y^j^ to represent the 

set {y®, i - 1, 2 , • • • , /:), in particular, y^' and y® denote yj and yj in the following examples for 
notational simplicity. |Y| = m, \ ■ \ denotes the cardinality of a set. 

We give a new index reduction technique for DAEs and define the weak differentiation index. 
With loss of generality, consider n polynomially DAEs with m dependent variables yj = y ft) with 
t a scalar independent variable, of the form 

f — f(t, the yj and derivatives of them) = 0, !</<«, 1 < j < m. (4) 


* where ^ denotes the transposition, which is the same way for the rest of this article. 


4 






It is the following equivalent form from the above notations, 

h 

fi = Cio + '^CikFik- ( 5 ) 

k=i 

where c,o, Cik are the coefficients that are the known forcing functions with t or the constants, 
P,vt = is a monomial in \y\,y 2 , • • ■ \ • • • ,y^m , ■ • • ,ym’'^] with exponent 

vector Oik and I, - \aik\ > 1. In particular, the highest degree of yj and its derivative y^’ denote 
dj and djr. in {/i,/ 2 , • • • ,/,), respectively. 

In order to compute the differential algebraic resultant in Section [H the outline of index 
reduction procedure is as follows: 

Phase 1 Initialization. 

f ^ 1 

(a) Collect every dependent variable yj and its derivative y.' for each fi, and then gather 
the set of dependent variables V. 

(b) Sort V for every yj and y ^' into a lexicographic order under assumption of ordering 

■ ■ ■ > y 2 > y^p^ > ■ ■ ■ > > yi, and Yj represents the set of yj and its derivative 

(c) Construct a matrix M - (niij), which is called variable pencil, defined for (|4]i by 

(1 the y j or its derivative appears in equation f, 

niij = < •' ( 6 ) 

10 if the variable does not occur, 

where row{M) and col{M) denote the number of rows and columns of M, respectively. 
Phase 2 Differentiation. 

(a) Determine the set of differential equations Fg if there exists niij - 1 for any derivative 
of Y, and the set of algebraic equations Fa if niij — 0 for all derivatives of Y in f, where 
\Fo\ = i, \Fa\ ^n- s. 

(b) Select the algebraic constraints fk{s +1 < k < n) from Fa to differentiate Vk such 
that ord(yj) < rj, which can be viewed as the homogeneous order. If it generates the 
new differential dependent variables, then it requires to augment the row and column of 
variable pencil to denote A1', update dependent variables set to V' and V'y The terminated 
condition of algebraic differentiation is as follows: 

n 

n+ |V'| - |V';| + 1. That is, row(M') = col(M') - (7) 

/:=.?+! 

where vi - V 2 - ■ ■ ■ - Vs - Q. 

(c) Select some low order differential equations fk{l < k < s) from Fg to differentiate 
Vk if (01 fails such that ord(fk) < maxJ^j rj with yj, and augment the row and column of 
variable pencil to denote A1", and update dependent variables set to V", V"y and the order 
r, to r'.. The terminated condition of differentiation is as follows: 

J 

n 

n + Yj k>k - |V"| - |V"y| + 1. That is, row(AI") = col{M”) - (8) 

k=i 
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Remark 3.1. We remark that the termination of index reduction procedure is required by the 
condition 0 or 0 because of the construction of a square elimination matrix. In particular, the 
procedure may have fully been degenerated into the problems ofODEs if fails in Phase 2(c ). 
We can obtain the new set of ODEs from Phase 2(b) and (c), then refer to the algorithm of U4l . 

Definition 3.1. The number of differentiations specified by index reduction procedure gives a 
formula for the weak differentiation index of system ofDAEs, denoted by d„ — max"^j Vk. Obvi¬ 
ously, if no differentiation of the original system is index zero, ODEs may have the weak differ¬ 
entiation index more than zero. 


Theorem 3.2. Let F = {/i,/ 2 , ■ • ■ ,/„} e Rfyuyi, ■ ■ ■ ,ym,y^l\ ■ ■ ■ ,y^l"\ ■ ■ ■ ,ym\ ■ ■ ■ ,ym"\ and 
the index reduction procedure satisfies the terminated condition 0 or 0- Then v — (ui, V 2 , - ■ ■ , Vn) 
can be computed correctly as specified. 

Proof. The initialization in Phase 1, we can easily get the set of dependent variables V and 
construct the variable pencil M - (mij) from F and initialize ui - U 2 - ■ ■ ■ - v„ - 0. According 
to Phase 2(a), we have 

{ /l,/2, ,fs fs+l,fs+2,- .,/« }■ 

ODEs (Fo) algebraic equations (F^) 

To compute u - {vi,V 2 , ■ ■ ■ , v„), two cases are considered: 

Case (a): only differentiate Fa - {/j+i, fs+ 2 , • • • , /«} to satisfy the condition 0 such that ordiyj) < 
rj based on the homogeneous order, which can repeat the differentiation to obtain the dif- 
ferentiation times {l'j+i,! 25 + 2 . ■ • ■ ,yn}- In the general setting, since |y^.^ | = rj, it is easy to 
get the terminated condition 0. In particular, since < rj for sparse case, it always 

generates the new differential dependent variables {yf'\y^ 2 ^\ ■ ■ ■ with kj < rj, and re¬ 

quires to update the set of dependent variables V' = V U U y^^^ U ■ ■ ■ U and F = 

F U {d/s+i, ■ • ■ ,d'''+‘/j+i, ■ • ■ ,<5/„, • ■ • ,d‘'''fn}. Consequently, we need to augment the row and 
column of variable pencil A1 to denote Af'. This concludes following: 

row(Ai') = col(M' \ \yj,y^j‘"^]) + 1- 

Case (b): following Case (a), if the condition 0 fails, it needs to obtain the condition ®. The 
problem can be transferred into the general nxm system of ODEs Eg - , • ■ • , 

In order to reduce the redundancy differentiation times, we only differentiate some low order 
ODEs from Eg such that ord(fk) < rj with yj, which are mjj = 0 in Af' for {yf‘\y 2 ^"\ • ■ • 
Eurthermore, it may need to differentiate some general ODEs from Eg to satisfy the condition 0, 
which can repeat the differentiation to obtain the differentiation times {vi,V 2 , ■ ■ ■ , v,,]. It always 
generates the new differential dependent variables {y ^, y ^^'^, • ■ • , y'tc '^} and requires to update the 

set of dependent variables V" = V'Uyj^'^Uy^^^'u- • -UyJ*'"^, F' = Fu{b/i, • • • ,6'’'f\, ■ ■ ■ ,5f„, • • • , 
and the order rj to r'.. Consequently, we need to augment the row and column of variable 
pencil AI' to denote At". This concludes following: 


row(At") = co/(AI" \ {yy.y^^'^) -i-1. 


□ 
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Remark 3.3. Yang et al n34l gave a formulation of differentiation times for differential elimi¬ 
nation ofODEs. However, their method may lead to some redundant differentiation times in the 
practical applications, such as the constrained mechanical systems. Here, we propose a variable 
pencil technique to analyze the differentiation times of DAEs. It is able to make differentiation 
times as few as possible for differential algebraic elimination. It is also suitable for the differen¬ 
tial elimination ofODEs with Dixon resultant formulation. 


We present a simple example to illustrate our index reduction procedure as follows: 


Example 3.1. This example is the linear, time-dependent index two DAE discussed in Gear Ml 
as follows: 


1 qt 


0 0 



1 +?7 pi\ 
dt J\y2/ 



( 10 ) 


where the dependent variables yi,y 2 with t a scalar independent variable, pi and p 2 are the 
known forcing functions oft, and q is a parameter. We can get the expanded form as follows: 


0 = /i = ji + (1 + q)y2 + qty2 - pi, 

0 = /2 = yi + qty2 - P2- 

We can initialize the original system (|77} as follows: 

(a) collect the set of dependent variable \ - {>'i,y 2 :>' 2 :yi}/ (b) sort the setY — [y\,y\,y 2 -,y 2 ]i 
(c) construct the variable pencil 



yi yi y 2 y 2 



h 

0 

1 

1 

1 


h 

1 

0 

1 

0 

Obviously, we can get the Fa and F^ 

with \Fa\ 

= \Fo 

1 = 

1 , ai 

homogeneous order as follows: 


yi 

y\ 

y2 

h 


/i 

ro 

1 

1 

1 

M ^ 

fi 

.1 

0 

1 

0 


Sf2 

0 

0 

1 

0 


Therefore, we have 


0 = /3 = 5/2 = ji + qy2 + qty2 - P2- 


( 12 ) 


For the differentiated equation 5/2 appended to the original system, the system of three equa¬ 
tions fi, f 2 and /s has four dependent variables y\,y\,y 2 , and y 2 . For eliminating the dependent 
variables (yi, ji) or {y 2 ,>' 2 ), the terminated condition of algebraic differentiation ^ holds. Con¬ 
sequently, we can get the differentiation times v — (0,1). 


Remark 3.4. We only differentiate the equation f 2 once, i.e., d„ — 1, and mix the algebraic 
equations and DDEs to deal with uniformly. However, the existing methods need to differentiate 
f 2 twice until no algebraic equations appear by substitution, that is, the problem is index two 
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4. Differential algebraic elimination 

In this section, the definition of differential algebraic elimination of DAEs is introduced by 
using the generalized Dixon resultant formulation. Based on the index reduction algorithm in 
Section 12 we also present an algorithm for computing the differential algebraic resultant. More¬ 
over, its basic properties are given. 


4.1. Definition of differential algebraic elimination 

The fundamental tool is based on the idea of algebraic Dixon resultant to create the differen¬ 
tial algebraic elimination. Firstly, we construct the differential algebraic cancellation matrix, and 
then compute the entries of differential algebraic elimination matrix, determinant of which con¬ 
tains the differential algebraic resultant as a factor. That is, DAEs can be treated as polynomial 
system, and yj and its derivatives can be viewed as parameters, the other dependent variables and 
their derivatives as the purely algebraic variables are eliminated simultaneously. Therefore, we 
can obtain the single ODE with yj and its derivatives to directly apply the numerical method. 


Let F = {/i, / 2 , •••,/„, (5/i(5/;,, , 6’^"fn] e /?[yi,yi, ■ • ■ ,y™,• • • ,yf'\ 


Tn 


Y = (yi, y 2 , • • • , y,„, y'/^ 


the system 


(/i =0,/2 = 0, •■•,/„ =0} 


(13) 


has solution if and only if the system F has solutions for v - {v\,V 2 , - ■ ■ , VnV- In order to 
define the differential algebraic elimination of (fT3l l it is necessary to find a weak differenti¬ 
ation index u for eliminating the yi,y2, • • • ,yj-i,yj+i, yj+i, ■ ■ ■ ,ym and their derivatives, such 
that /i, • ■ • ,fn,6f\, ■ ■ ■ ,6'^'fi, ■ ■ ■ ,6'^"f„ are {v\ -i- 1) -i- (u 2 H- 1) H- ■ ■ • H- (l»„ H- 1) polynomials in 
vi + V 2 + ■ ■ ■ + v„ + n - 1 variables. 

By following the definition of Dixon resultant we have 

Definition 4.1. Let fi be a differential polynomial in R[yi,y 2 , ■■ ■ ,ym,y^i \ ''' ?y^('^\''' jy»!\''' jym”^]> 
N — (l'i-I-1)h-(l'2-I-1)h-- ■ ■-l-(i/„H-l)— 1 as mentioned above. The differential algebraic cancellation 
matrix DCw of¥ with yj and its derivatives is the (A H- 1) X (A + 1) matrix as follows: 



fi(yi,y 2 ,- 



in 

• .ym 

^) • 

/iv+i(yi,y2, • 

• AmOl'V 


9 jm 


fiiyuyi,- 



in 

^) • 

•• /iv+i(yi,y2, • 



■ ,yy 

DCw - 

fiiyuyi,- 



in 

• ,ym 

■^) • 

/iv+i(yi,y2, • 



■ yy 


.fi(yi,y 2 , ■ 



-in 
• ,ym 

'^) • 

/w+i(yi-y2, • 

■ ,ym,y^i\- 

■yyr 

■ yy 


where {y^.y^^ • • ■ tisparameters do not replace by {yj,y^j \ • • • bi /i(l < / < A H- 1). 

The differential algebraic elimination polynomial of¥ is denoted by DO^ e /J[Y, Y], 


DO^ 


det(!DCw) 




(14) 
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the row vector of differential algebraic elimination derived polynomials o/F is denoted by DP^, 
and differential algebraic elimination matrix o/F is denoted by DAf as follows. 


D6f - 




di-\ 

yt 


j 


■DAf 




3^1 


n“i 

i*j 




where the rows and columns of DA^ are indexed ordering by y\ 

■■■ > Ji’’ > ym > ■ ■ ■ > y\, fw"'’ > 




> yi? > ■ ■ ■ > > 


(15) 

> > ym > ■■■ > yi. 


respectively. The coefficient matrix DA^ is also a square matrix, determinant of which is called 
differential algebraic resultant, denoted by DARes{yj,y ). 


Here, we can write the DAf in the following block structure notation; 


DAk - 


£>0,0 

£>1,0 


Dd, 


1,0 


£>0,1 

£>1,1 

£>i/,-i, 


£>0,OTi-1 

£>1,^1-! 


Ddy 


l,lVrf,-lJ 


(16) 


where each block Djj is of size {N - rj - 1)! Yiti n);,=i x (N - rj - 2)! Yiti n^;=i didiin- 

i*i ' i*i 

As the increasingly large scale system of DAEs, we can make use of its structure and block 
triangularization to decompose a problem into subproblems by permuting the rows and columns 
of a rectangular or square, unsymmetric matrix. For more details refer to 12211 . 

Following the properties of Dixon resultant we prove easily. 


Theorem 4.1. The differential algebraic elimination matrix DA^ is of size N\ Y\T=i didif^. x 

i*i 

N\ nr=i n:/ .j djdif^. at most, and the complexity of computing the entries ofDAr is 0{d^{N\ Y\T =2 
i*i ' i*j 

didi/j.ff) in the worst case, where di and di^, are mentioned above. 


Theorem 4.2. Differential algebraic resultant can be expressed as a linear combination of en¬ 
larged system of equations F with yj. 


II Uj 

DARes(yj,y]^^^) = Z 
(=1 ^,=0 


(17) 


where is a polynomial with respect to Y and can be deduced from DP^. Moreover, if 
DARes{yj,yj‘ ) is a reducible differential equation, it can also be proved that the extraneous 
factors mentioned above may include three parts which are taken from DPf, DA^ and the result¬ 
ing resultant expression by substituting DP^, respectively. 
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Remark 4.3. From Theorem 14.21 we can remove the extraneous factors from differential alge¬ 
braic resultant when the existing greatest common divisors in each row or column of differential 
algebraic elimination matrix. That is, the extraneous factors are the greatest common divisors in 
the algebraic cofactors ofDAf. 


Theorem 4.4. Differential algebraic resultant is equal to zero that is a necessary condition for 
the existence of a common solution of system ofDAEs. 

Corollary 4.5. Let yi,y 2 , - ■ ■ ,ym be solutions of the system ofDAEs (El- Then yj satisfies the 
DARes(yj,y^l\ 


4.2. Algorithm 

In this subsection, we have the following procedure for differential algebraic elimination. 


Input; DAEs system F — {f\ = 0,/2 = 0, • • • ,/„ = 0}, and dependent variables Y \ {yj,y^j^^]. 
Output: a polynomial ODE only contains yj and its derivatives. 


rfi 

Step 1: Count the number of DAEs and Y \ {yjO/ !’ denote n and m respectively, if n is equal 
to m plus 1, then goto Step 3. 


Step 2: Call index reduction algorithm in Section|3]by taking the f-derivative of f, update n and 

1 

m such that n = m + 1 by the enlarged system of equations F and new Y \ 1 ’ the 

collections are as follows. 


fudfu 

/2,b/2, 

■■,6^'f, 

■ ■ , 

■ =0, 

Y «- • 



Jn,6fn, 




Am. ylnK 

■■ 


\ } 


( 18 ) 


Step 3: Construct the differential algebraic cancellation matrix !DCw, obtain the entries of differ¬ 
ential algebraic elimination matrix DAr, remove the greatest common divisors from each 
row or column of DA^, and then compute its determinant DARes(yj,y .‘). 

Step 4: Return DARes{yj,y^^d^). 

Theorem 4.6. The above algorithm works correctly as specified and its complexity mainly con¬ 
tains index reduction algorithm and the computation of differential algebraic resultant. 


Proof. Correctness of the algorithm follows from the Dixon elimination method. Regarding 
the dependent variables {yj, y. ‘ ) as parameters and the other ones as algebraic variables, we can 
treat the enlarged system of equations F as an algebraic system. As shown in ll34ll . a necessary 
condition for the existence of a common solution of algebraic differential equations is that the 
differential resultant is equal to zero. We can easily get the Theorem l4.4l 

Erom the description of algorithm, we observe that there are two major steps on time com¬ 
plexity. In Step 2, we can obtain the differentiation times YIk=i ^k- The problem is solved by 
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homogeneous order rule, one that makes differentiation time as few as possible for reducing the 
number of enlarged system of equations F. If there exists the v - {vi,V 2 , - ■ • , VnY, which can be 
done in polynomial time. In Step 3, the complexity includes three parts: (a), to obtain the entries 
of differential algebraic elimination matrix DA^ in Theorem l4.ll suppose = d, ri - r, 

for each i — 1 , 2 , • • • , m, it needs at most 

m r, m r,- 

Y] n n n 

i-2 /i, = l 1=1 /i,-=l 

which is the single exponential complexity; (b). calculate the greatest common divisors for each 
row or column of DA^ in the polynomial time; (c). compute its determinant DARes(yj,y d ) 
in polynomial time. Therefore, if there exists the differential algebraic resultant with single 
dependent variable and its derivative, we can transform the system of DAEs to its equivalent 
ODEs in the single exponential complexity. □ 

Remark 4.7. From Lemma [Ol awt/ Theorem \4.6\ our algorithm is not a complete method. How¬ 
ever, our algorithm is really effective and practical technique on numerous problems. It is well 
known that Dixon resultant elimination is that it can do one-step elimination of a block of un¬ 
knowns from a system of polynomial equations like Macaulay’s. Moreover, the size of Dixon 
matrix is much smaller than the size of Macaulay matrix. Though the entries of Dixon matrix 
are complicated in contrast to the entries in Macaulay matrix, the entries of which are either 
0 or coefficients of the monomials in the polynomial systems. Fortunately, we can easily apply 
the extended fast algorithm for constructing the Dixon matrix SMi- In particular, for a fixed 
number n of variables of a polynomial system, the construction cost of Dixon matrix is at most 
0{mvol(fff) arithmetic operations where mvol{7^) is the n-fold mixed volume. As shown 
in & our algorithm is also appropriate to the system of ODEs. In many practical applications 
of DAEs, we can easily see that di and di^. are very low degrees in 

Example 4.1. Continue from Example \3.1\ we can construct the differential algebraic elimina¬ 
tion matrix with yi and y\ as follows: 

[qtpi -qtp 2 +yi . 

We can also eliminate yi and ji by the same method simultaneously. Therefore, we get the 
following differential algebraic resultants: 

DARes(yi,yi) = yi - /72 + r]t(Pi “ P 2 ), 

DARes(y 2 ,y 2 ) = yi - Pi + Pi- 

These are the same as the results obtained in ini/- 

5 . Examples 

In this section, we present some small examples in detail and compare the matrices size of 
differential resultants of two models with other methods. Examples 15 .1 1 and 15.41 illuminate how 
to deal with the nonlinear and high index DAEs of constrained mechanical system. Example l5.2l 
uses a simple example to test our algorithm for the nonlinear and non-square system of DAEs. 
Example 15. 3 1 is a practical application for the linear electrical network problem. 
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5 . 1 . Some examples in detail 

Example 5.1. Consider the nonlinear DAEs for the simulation of the dynamics of multibody 
systems, which is a major application area. Here, we show the simple pendulum to illustrate 
many of the principles. The DAEs can be written 


0 = /i =yi 
0 ^ fi^yi+yi^- g, 




( 19 ) 


where g > 0 ,L > 0 are constants. From its variable pencil, labeled by equations and 
variables, is 




fi 

h 


yi 

1 

0 

1 


yi yi 
1 0 
0 1 
H] 1 


h 

0 

1 

"o 


A 

1 

1 

0 


Obviously, we can get the Fa and F 

0 with |Ea| 


= 2, 

and differentiate fj, based 

homogeneous order as follows: 









y\ 

ji 

y\ 

y 2 

n 

h 

A 

/i 

1 

0 

1 

0 

0 

0 

1 


/2 

0 

0 

0 

1 

1 

1 

1 


M h 

1 

0 

0 

1 

0 

0 

0 


Sfi 

1 

0 

0 

1 

0 

0 

0 


1 

1 

0 

1 

1 

0 

0 


Therefore, we have 

0 = /4 = 5/3 = 2 yi + 2 y 2 y 2 , 

0 = /5 = 5^/3 = 2 yiyi + 2y2% + 2 y\ + 2 y\. 

For the differentiated equations 5/3 and d^fj, appended to the original system, the system of five 
equations /i,/2,/3,/4 and f$ has seven dependent variables y\,yi,y\,y2,y2,y2, and A. For elim¬ 
inating the dependent variables {yi.ji.yi) or \y2,y2,y2], ths terminated condition of algebraic 
differentiation satisfies Consequently, we can get the differentiation times v — ( 0 , 0 , 2 ), that 
is, d„ - 2 . 

We can construct the differential algebraic elimination matrix with , y\ and yi as follows: 



y\y\ 

0 

0 

1 

0 

0 

0 


L^-y\ 0 0 

0 -ylh-yifi yiyi 

0 L^-y] 

0 0 0 

-1 -gyi 0 

0 -yi 0 

0 0 0 


0 0 0 

y]-L^ 0 0 

0 0 0 

0 0 

0 F^-y] -yjyi - yiyj 

1 0 -gyi 

0 -1 -yi 
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Obviously, we can also eliminate yijji and yi by the same method simultaneously. Therefore, 
we get the following differential algebraic resultants: 


DARes(yuy\,yi) = + y\L'')y\ + {-IL'^y] + 2L^yy)y\yy + 3g^L"'y\ + + g^y\ 

-3g^Ly,-L^gy„ 

DARes(y 2 ,y 2 ,h) = - gy\ + 2 gLy - gL^. 


The remaining dependent variable A is determined by yi and y^- 


Example 5.2. The example is the nonlinear, non-square system of DAEs discussed in iz3/ as 
follows: 


Cio 0 
C20 0 

C30 C 31 


0 

C 22 

0 


Cl3 

0 

0 


1 'I 

ym 

ym 


0 

0 

0 


( 21 ) 


where the dependent variables y\ and y^, Cjjii — 1,2, 3,7 = 0,1,2,3) are the known forcing 
functions oft. We can get the expanded form as follows: 


0 = /i = cio + Ci3yiy2, 
0 = /2 = C 20 + C22yiy2, 
0 = /3 = C30 + C 3 iyiy 2 - 


( 22 ) 


We can initialize the original system ( 1221 ) as follows: 

(a) collect the set of differential variable V = (yi,y 2 ,>' 2 .>'i}/ (b) sort the set V = {yi,yi,y 2 ,y 2 }/ 
(c) construct the variable pencil 




fi 

f2 

h 


yi 

0 

0 

_i 


ji y2 

1 0 

1 1 

1 )\ 1 


y2 

1 

0 

0 


Obviously, we can get the Fa and Eg with |Ea| = 1 and |Eo| = 2, and the system of three equa¬ 
tions fi, f 2 and fi, has four dependent variables yi,yi,y 2 ‘^nd fi. For eliminating the dependent 
variables {yi,yi) or {y 2 ,y 2 }, the terminated condition of algebraic differentiation satisfies 0 - 
Consequently, we can get the differentiation times v — (0,0), that is, d„ — 0. 

Here, we can construct the differential algebraic elimination matrix with yi and yi as follows: 

[c20C3iyi - C22C30yi]j^j ■ 

We can eliminate yi and yi by the same method simultaneously. Therefore, we get the following 
differential algebraic resultants: 

DARes{yi,yi) = C 2 oC 3 iyi -C 22 C 3 oyi, 

DARes(y 2 ,fi) = -Cl 0 C 22 y 2 +C 20 Cl 3 y 2 - 
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Example 5.3. Consider a practical linear electrical network example, differential algebraic 
equations of index 1 may have an arbitrarily high structural index from as follows: 


0 = /i ^yi + h+y\- a(0, 
0^ fi^yi + h+ yi- b(t), 
0 = /3 = j4 + js + ^3 - c{t), 
0 = /4 = ^4 + js + y4 - d{t), 
0 = fs^ys- eit), 


(23) 


where a(t), b(t), c(t), d(t) and e(t) are the known forcing functions oft. It is clear thatys is known, 
i.e., ys — e{t). We can get the variable pencil as follows: 


fi 

M => /3 

/4 

h 


yi yi yi ys 

10 10 
0 110 
0 0 0 1 

0 0 0 0 

0 0 0 0 


^3 y4 ^4 ys 

10 0 0 

10 0 0 

0 0 10 
0 110 
0 0 0 1 


ys 

0 

0 

1 

1 


Obviously, we can get the Fa and Fa 

with \Fa 


lAFo 

1 = 4, 

and differentiate fs based on the 

homogeneous order as follows: 











yi 

yi 

h 

ya 

h 

y 4 

y 4 

ys 

ys 

fi 

rl 

0 

1 

0 

1 

0 

0 

0 

0 


h 

0 

1 

1 

0 

1 

0 

0 

0 

0 


M! ^ 

0 

0 

0 

1 

0 

0 

1 

0 

1 


14 

0 

0 

0 

0 

0 

1 

1 

0 

1 


fs 

0 

0 

0 

0 

0 

0 

0 

1 

0 


Sfs 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Therefore, we have 











0 


II 

II 

ys - 

e(t). 




(24) 


For the differentiated equations 5fs appended to the original system, the system of six equations 
/i>/ 2 ./ 3 ./ 4,/5 and ff has seven dependent variables yi,y2,y2,y3,y3,y4 andy^. For eliminating 
the dependent variables {y 4 , ^ 4 ), ty 3 ,>' 3 } or {y 2 ,>' 2 !. the terminated condition of algebraic differ¬ 
entiation satisfies 0 - Consequently, we can get the differentiation times v — (0,0,0,0,1). It 
only needs to differentiate fs once, that is, d„ — 1. 

Remark 5.1. We easily compute the differential times (0,0, 1, 1,2) of five equations ff, f 2 , f 3 , ff, fs 
in sequence by the 'L-method that is, the structural index is 3. As for an increasingly large 
dimensions, the 1,-method may perform an arbitrarily high differentiation times. However, our 
weak differentiation index is the same as the differentiation index. In general, it is suitable for 
the linear DAEs as follows, 

AY + BY = S, (25) 
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where S is the vector of m sufficiently smooth forcing functions oft, Y is as mentioned above, A 
and B are m X m matrices, such as linear DAEs ( 1251) with m — 2k + 1, zero vector §, and the 
identity matrix B, 


1 1 

1 1 

■ 1 1 
1 1 


1 1 
1 1 


such that A solely consists ofk blocks of form ( [ [ ), the lower left element of each being on the 
main diag onal of A. This is the same result that the index is 1 by using the Kronecker canonical 
form \3l\l . However, structural index algorithm Ml needs to differentiate the last equation k 
times, that is, the structural index is k + \. In Ml. Pantelides’ algorithm needs to perform k + \ 
iterations before termination. Therefore, it leads to a large number of redundant differentiation 
times. 

Example 5.4. We present a double pendulum model to demonstrate our index reduction tech¬ 
nique in the dynamical systems. It is modeled by the motion in Cartesian coordinates, see Figure 
I. 


( 0 = 0 ) 



Figure 1: Double Pendulum 
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We can derive the governing DAEs using Newton’s second law of motion as follows, 

0 = /i = mixi - - ^(X 2 - xi), 

0 = /2 = miyi - ^yi - ^(j 2 - yi) - mg, 

0 = /3 = m2X2 - ^{X2 - Xi), , ^26) 

0 = /4 = m 2 y 2 - ^(yi -yi)- m 2 g, 
o^fs^^x+y\-ii 

0 = fe ^ (X 2 -Xif+ (y 2 -yif - ll, 

where g > 0,li > 0 ,l 2 > 0,mi > 0 and m 2 > 0 are constants, the dependent variables 
, x, 2 , y \, y 2 with t a scalar independent variable, Ai, A 2 are the known forcing functions oft. 
From m, we can construct its variable pencil M, and easily to get the needed to differenti¬ 
ate fs and ff twice, respectively. Then we obtain the new variable pencil M' to determine the 
termination of differentiation, which holds the condition 0- Finally, we can eliminate the depen¬ 
dent variables {xi, xi, xi}, {x 2 ,X 2 , X 2 }, {yi,yi,yi}, or \y 2 ,y 2 ,y 2 ], where Ai,A 2 can be determined 
by xi,X 2 ,yi and y 2 . Therefore, we can get the differentiation times v - (0,0,0,0,2,2), that is, 
d„ - 2. 

5.2. Some comparisons 

In this subsection, we also apply our algorithm to the system of ODEs and compare the 
matrix size of differential resultant with other methods. The algebraic manipulation of differential 
equations has gained importance in last years. 

Example 5.5. Consider two nonlinear generic ordinary differential polynomials with order one 
and degree two from /Hz as follows: 

0 = /i = j? + aiyiyi + a2y] + a^yi + a4yi + as, 

0 - /z = ji + ^lyiji + b2y\ + b^yi + b^yi + bs, 
where a,-, h, are differential constants, i.e., da, — dbj = 0(/ = 1,2, • ■ • , 5). 

Example 5.6. Consider the simplified version of a predator-prey model from id/ as follows: 

0 = /i = a2yi + (ai + a4yi)y2 +y2 + (as + aeyi)yl + asyl, 

0 = /2 = ji + (hi + b3yi)y2 + (ha + hsyOyj + h4y2, 
where a,, bj{i = 1,2, ■ • ■ ,6,j- 1,2, • • ■ ,5) are the known forcing functions of t. 




Example 

Matrix size 

ZYG [3^ 

Rueda [3ff] 

Our algorithm 

15.51 

36 X 36 

* 

9x9 

EE 

* 

13 X 13 

5x5 


Table 1: Matrix size for computing dilferential resultant 
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Table 1 gives a comparison of the matrix size of differential resultant in Examples 15.51 and 
15.61 where represents that the computation is not compared. From the Table 1, we have the 
following observations: 

In two examples above, the matrix size of differential resultant via our algorithm is much 
smaller than two other methods. The smaller matrix leads to reduce more time for computing 
its symbolic matrix. It is consistent with the generalized Dixon resultant formulation. ZYG 
is based on the idea of algebraic sparse resultant and Macaulay resultant for a class of the 
special ordinary differential polynomials. Rueda S presents the differential elimination by 
differential specialization of Sylvester style matrices to focus on the sparsity with respect to the 
order of derivation. In practice, Dixon’s method is the most efhcient technique to simultaneously 
eliminate several variables from a system of nonhomogeneous polynomial equations. 

6 . Conclusions 

In this paper, we propose a new index reduction for high index DAEs and establish a relation¬ 
ship between the generalized Dixon resultant formulation and system of DAEs solving, which is 
defined as differential algebraic elimination. A significant problem in the differential algebraic 
elimination is to create methods to control the growth of differentiations. Our method can be 
applied to the mixed algebraic equations and differential equations to deal with simultaneously, 
and given a variable pencil technique to determine the termination of differentiation. From the 
algebraic geometry, it can be considered as the index reduction via symbolic computation. 

Our method can be also suitable for the system of ODEs and the high index nonlinear non¬ 
square system of DAEs, i.e., the number of dependent variable is not equal to the number of 
equations. The weak differentiation index is defined to unify the formulation of differentia¬ 
tion times for differential elimination of ODEs and differential algebraic elimination of DAEs. 
Moreover, a heuristics method is given for sifting the extraneous factors in differential algebraic 
resultants to remedy the drawback of factoring large polynomial system. Parallel computation 
can be used to speed up the computation of differential algebraic resultant of each dependent 
variable. 

However, the disadvantages of our method contain its limitation to polynomial coefficients 
and incomplete method because of the generalized Dixon elimination. Usually, for many prac¬ 
tical relevant applications, the large scale system of ODEs/DAEs is also a challenge problem by 
the purely symbolic method; for instance, the full robot in the Modelica context 191] has before 
symbolic simplification about 2391 equations and 254 dependent variables, which are reduced 
to 743 equations and 36 states that require a lot of index reduction going on. An obvious future 
work, is to attempt the block triangularization and sparsity considerations in constructing the 
differential algebraic elimination matrices. The sparseness is reflected in the quantity /,■ of P„t 
in Section |3] Furthermore, symbolic-numeric differential algebraic elimination method is a very 
interesting work in the numerical algebraic geometry. 
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